	PRO ENERGY_MASSIVE,dir,Nproc,ftime,PS=PS
; Calculate energy from massive files
; ftime is the final frame number (scalar)
leer_input,dir,N,nu,eta,eps,delta,dt,nstep,nstg,nstdat,nstk,nstp
x=2.*!pi*dindgen(N)/N & y=x
dx=x(1)

tiempo=dblarr(ftime+1)
Kp = tiempo ; kinetic energy proton
Ke = tiempo ; kinetic energy electron
EB = tiempo ; magnetic energy

; For each frame
for i=0,ftime do begin
time ='000'+string(i,format='(i1)')
if i gt 9 then time ='00'+string(i,format='(i2)')
leer,dir,time,Nproc,a,f,bz,uhz,jhz,whz,t
comcampos,f,uhx,uhy,x,y ; Compute hydrodynamic velocity components
comcampos,a,bx,by,x,y ; Compute magnetic field components
comcampos,bz,jhx,jhy,x,y ; Compute hydrodynamics electric current components
tiempo(i) = t
Kp(i) =  (dx^2./2.)*total(uhx^2.+uhy^2.+uhz^2. + ((eps*delta)^2./(1+delta)^2.)*(jhx^2.+jhy^2.+jhz^2.) + 2.*delta*(uhx*jhx+uhy*jhy+uhz*jhz)/sqrt(1+delta))
Ke(i) =  delta*(dx^2./2.)*total(uhx^2.+uhy^2.+uhz^2. + (eps^2./(1+delta)^2.)*(jhx^2.+jhy^2.+jhz^2.) - 2.*eps*(uhx*jhx+uhy*jhy+uhz*jhz)/sqrt(1+delta))
EB(i) =  (dx^2./2.)*total(bx^2.+by^2.+bz^2.)
endfor
Et = Kp + Ke + EB
!p.multi = [0,1,1]
if not keyword_set(ps) then window,0,xs=600,ys=600
if keyword_set(ps) then begin
set_plot,'PS'
device,filename=dir+'energies_massive.eps',bits_per_pixel=8,/encapsulated,xs=15,ys=12
endif
plot,tiempo,Eb-Eb(0),linestyle=3,xr=[0,max(tiempo)],yr=[-0.005,0.005],tit='Total (solid) Magnetic (.-) Proton (-) Electron (.)',ytit='Energy',charsize=2.0
oplot,tiempo,Ke-Ke(0),linestyle=1
oplot,tiempo,Kp-Kp(0),linestyle=2
oplot,tiempo,Et-Et(0),linestyle=0
if keyword_set(ps) then begin
device,/close
set_plot,'X'
endif
end
